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A surface of a strong topological insulator (STI) is characterized by an odd number of linearly 
dispersing gapless electronic surface states. It is well known that such a surface cannot be described 
by an effective two-dimensional lattice model (without breaking the time-reversal symmetry), which 
often hampers theoretical efforts to quantitatively understand some of the properties of such surfaces, 
including the effect of strong disorder, interactions and various symmetry-breaking instabilities. Here 
we formulate a lattice model that can be used to describe a pair of STI surfaces and has an odd 
number of Dirac fermion states with wavefunctions localized on each surface. The Hamiltonian 
consists of two planar tight-binding models with spin-orbit coupling, representing the two surfaces, 
weakly coupled by terms that remove the extra Dirac points from the low-energy spectrum. We 
illustrate the utility of this model by studying the magnetic and exciton instabilities of the STI 
surface state driven by short-range repulsive interactions and show that this leads to results that 
are consistent with calculations based on the continuum model as well as three-dimensional lattice 
models. We expect the model introduced in this work to be widely applicable to studies of surface 
phenomena in STIs. 



I. INTRODUCTION 

The physical feature that makes topological insu- 
lators special is their surface statesP^ The surface 
states are topologically protected, have been well doc- 
umented experimentally^ 4 -"^ and are predicted to exhibit 
a wide range of interesting phenomena when subjected 
to various perturbations?"^ Also, they harbor a poten- 
tial for future practical applic ations in spintronicsj 16 * 17 ! 
low-dissipati on ele ctronics^ 1 ^ and topological quantum 
computation! 2212 ^ 

The most interesting in this regard are the three di- 
mensional strong topological insulators whose surface 
hosts an odd number of gapless fermionic states, while 
their bulk remains fully gapped. Of these, materials in 
the Bi2Se3 family represent a canonical example with a 
single Dirac state located at the center of the surface Bril- 
louin zone. At low energies, this state can be described 
by a continuum Dirac Hamiltonian, 

hu = v(k x s y - k y s x ), (1) 

where s = (s x , s y ,s z ) are Pauli matrices in the spin space 
and v is the characteristic Fermi velocity. The spectrum 
is ek = ±u|k| and can be gapped only by adding a term 
proportional to s z to the Hamiltonian, thus breaking the 
time-reversal symmetry (T). 

Now suppose we wish to describe this surface state by 
a 2D lattice model such that the above Dirac spectrum 
would emerge at low energies. Discrctizing Hamiltonian 
e.g. on a simple square lattice we obtain a candidate 
Hamiltonian 

/ijf u = v(s y sin k x — s x sin k y ) , (2) 

with the spectrum e^ u = ±« y sin 2 k x + sin 2 k y . We ob- 
serve that this spectrum indeed has a Dirac point at 
k = (0,0) but additional 3 Dirac points at (0,7r), (it, 0) 



and (7T, 7r) have been introduced by the discretization pro- 
cedure. While it is possible to remove these 'extra' Dirac 
points, e.g. by adding a term s z (2 — cosk x — cosk y ) to 
the Hamiltonian ^ , the price one pays is a Hamiltonian 
with broken T which therefore cannot faithfully describe 
the physics of the T-invariant state on the STI surface. 

The above construction illustrates the constraints im- 
posed by the well-known Nielsen-Ninomyia theorerrPH, 
which states that it is impossible, as a matter of princi- 
ple, to construct a T-invariant lattice Hamiltonian with 
an odd number of Dirac fermions in the low-energy spec- 
trum. The theorem assures us that the surface states of 
a STI do not have a faithful description in terms of a 
two-dimensional lattice model. To describe such surface 
states one either must use the continuum model Eq. ([I]), 
or, if a lattice description is required, one must consider 
a 3D lattice model with open boundary conditions along 
at least one spatial direction. 

The above conclusion imposes a number of restrictions 
on the theoretical studies of the phenomena associated 
with the STI surface states. While many qualitative prop- 
erties follow directly from the continuum effective theory 
the latter often fails when quantitatively accurate predic- 
tions are required. This is because quantities computed 
from the Hamiltonian typically depend strongly on 
the ultraviolet cutoff that must be imposed in order to 
regularize divergences. For instance, as explained in more 
detail below, the critical coupling for interaction-driven 
magnetic and excitonic instability in the Hamiltonian 0) 
depends linearly on the ultraviolet cutoff A, which leads 
to considerable uncertainty since the latter is not well 
known. Also, predictions that are strongly cutoff depen- 
dent must be viewed as inherently unreliable. To obtain a 
quantitatively reliable prediction one must resort to a 3D 
lattice model in a geometry with surfaces (e.g. a slab). In 
this case, calculations are typically to be performed nu- 
merically because of the low symmetry of the problem. 
Such calculations are often difficult, essentially because 
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the vast majority of the computational effort is spent on 
the bulk states which do not show any interesting physics 
but nevertheless must be included in order to capture the 
correct surface physics. 

In this paper we show that it is possible to construct 
a simple 2D lattice model that faithfully describes the 
physics of a pair of STI surfaces, such as those termi- 
nating a slab. The idea is based on a simple observa- 
tion that a pair of STI surfaces will have an even num- 
ber of Dirac points and therefore the Nielsen-Ninomyia 
theorem-- does not prevent us from constructing a rele- 
vant 2D lattice model. Our basic model-building strat- 
egy is straightforward: we describe each surface by a 
2D Hamiltonian very similar to Eq. ([2| and then con- 
struct a T-invariant coupling between the two surfaces 
that preserves the gapless Dirac points at k = (0, 0) but 
hybridizes, and therefore gaps out, the remaining Dirac 
points. This way, we end up with a 2D lattice Hamilto- 
nian describing a pair of parallel STI surfaces each con- 
taining a single gapless Dirac point. The important fea- 
ture of the model is that the low-energy eigenstates are 
localized on a single layer with a negligible overlap be- 
tween layers, and therefore faithfully capture the physics 
of the STI surfaces. 

To ascertain the validity of the proposed model we con- 
duct a number of tests. We show that applying vari- 
ous perturbations (such as inducing a magnetic gap at 
the Dirac point) to a single layer only affects the low- 
energy physics of that layer and leaves the other layer 
unchanged, exactly as one would expect in a physical 
STI sample. We also calculate the spectral function of 
the surface state from a full 3D model and demonstrate 
that, with a judicious choice of parameters, the 2D model 
introduced here gives the same result in the low-energy 
sector. 

Finally, as a potentially useful application of our new 
model we study the mean-field phase diagram of a thin- 
film STI in the presence of a short-ranged repulsive in- 
teraction (both intra- and inter-layer) and under exter- 
nal bias. This configuration exhibits two leading insta- 
bilities: a magnetic gap associated with spontaneous T 
breaking and an excitonic gap, driven by condensation 
of electron-hole pairs residing in different surfaces. Both 
the magnetic and the excitonic phase exhibit interesting 
physical properties. The magnetic phase has been pre- 
dicted to show the half- integer quantum Hall effect prim- 
age mag netic monopole effect™ and inverse spin-galvanic 
effect to name just a few, while the exciton condensate 
(EC) phase is expected to show dissipationless transport 
of electron-hole pairs, interesting thermo-electric prop- 
erties, as well as fractionally charged vortices!- 18 ! 24 " ^ In 
addition, the interplay between the magnetic and EC or- 
ders has been argued to exhibit a number of interesting 
phenomena, such as the anyon exchange statistics for cer- 
tain topological defectsp3 Similarly, interfacing such EC 
with a superconductor is predicted to result in interesting 
Majorana edge modesPS Using the proposed model we 
map out the corresponding phase diagram and study its 



evolution in response to the changes in parameters that 
can be tuned in a laboratory setting. The understanding 
that we thus gain will aid experimental searches for these 
interesting phases of quantum matter. 



II. LATTICE MODEL FOR THE SURFACE 

The basic strategy for our model design has been al- 
ready outlined in Sec. I. Here we provide the necessary 
technical details. To summarize our main objective we 
wish to construct a 2D lattice model whose low-energy 
properties will coincide with those of a surface of a re- 
alistic STI as modeled, e.g. by the standard 3D lattice 
model. To this end it is useful to visualize a STI in a slab 
geometry as composed of individual layers (e.g. the quin- 
tuple layers in Bi 2 Se3). We can now imagine integrating 
out the electronic degrees of freedom associated with the 
bulk layers, which are gapped, and retaining only the de- 
grees of freedom residing in the outermost surface layers. 
Such a procedure will leave us with a lattice model for 
the two surface layers, coupled through the bulk degrees 
of freedom that have been integrated out. 

Implementing this procedure in a lattice model turns 
out to be quite cumbersome and we defer this to a later 
section. Here, we proceed to construct the 2D model us- 
ing a simple heuristic procedure based on the symmetries 
and general considerations. 



A. Model I 

Our basic building block will be the effective 
HamiltoniarJ2122lfor a single quintuple layer of Bi 2 Se3 reg- 
ularized on a simple square lattice, 




where 

/ik = 2A(sj, sin k x — s x sin k y ), (4) 

and 

Mk = e — 2i(cosfc x + cos k y ). (5) 

The 2x2 matrix structure in Ho reflects the two orbital 
degrees of freedom per quintuple layer that are neces- 
sary to capture the physics of the topological phase in 
the Bi2Se3 family of materials. In the following we take 
A = 1 and measure all other energy scales in units of A. 
In the bulk the model parameters e and t have certain 
fixed values that are material specific. At the surface 
of a STI, however, we expect these parameters to be ef- 
fectively renormalized so that gapless states can emerge. 
Considering the spectrum of H it quickly becomes clear 
that the condition for the gapless states to occur at the 
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FIG. 1: Band structure of the effective surface Hamiltonian for model I (top row) and model II (bottom row). Panel (a) shows 
the spectrum of Hi (top) and Hn (bottom) in the first Brillouin zone along the path connecting points of high symmetry 
(-7T, tt) — > (0, 0) — > (0, 7r) — > (n, 7r). Panel (b) displays the effect of biasing the layers through the inclusion of SHv with V = 0.2 
and (c) shows the effect non-zero exchange coupling m = 0.2 in one of the surface layers. The parameters used are t = 0.5 and 
e = 2.0. 



r point [i.e. at k 
thus introduce 



(0, 0)] is e = it. For future use we 



Mk = 2t(2 — cos k x — cos k y ) 



(6) 



to denote this fine-tuned interlayer coupling. With this 
choice of parameters Mk=r vanishes and Hq will have 
two copies of a gapless Dirac fermion at this point, one 
associated with each orbital. If Hq is to describe a sin- 
gle STI surface then one of these must be removed. As 
discussed previously, the only way to accomplish this in 
a 7~-invariant manner is to introduce coupling to the de- 
grees of freedom associated with the other surface. This 
consideration inspires the following 8x8 Hamiltonian for 
the pair of surfaces, 
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(7) 



where f^k represents the aforementioned coupling. The 
Hamiltonian Hj, which we henceforth call 'model F, acts 
on an 8-component wavefunction where ^ i rep- 

resent the orbital and spin degrees of freedom associated 
with surface I — 1,2. There is a considerable freedom in 
selecting the form of the coupling between the surfaces, 
the only significant constraint being that f2 k 7^ near 
the r point, so that the redundant Dirac point is indeed 



removed from the spectrum. Thus, f^k = const is one 
possible choice, although we find that taking 
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[e + 2<(cos k x + cos k y 



(8) 



leads to better results when the overall spectrum is com- 
pared to the exact surface spectrum in a 3D STI model 
discussed below. 

The spectrum of Hamiltonian ^ consists of four dou- 
bly degenerate branches given by 



1/2 



(9) 



with e k = 4(sin 2 k x + sin 2 k y ). It is easy to see that near 
the r point, where Mk vanishes, the spectrum consists of 
four gapless branches, comprising two Dirac points, and 
four gapped branches. The full spectrum is displayed 
in Fig. [T^,. It is also easy to see that the wavefunctions 
associated with the two Dirac Fermions are localized on 
different surfaces: near the T point, where Mk vanishes, 
the two outermost diagonal elements in the Hamiltonian 
([7]) become decoupled from the rest of the system and 
thus the corresponding eigenstates can be chosen to have 
support in a single surface. 

It must be noted that slightly away from the T point 
Mk becomes non-zero and this leads to a weak mixing of 
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the two surface states. The strength of this mixing can 
be ascertained from Eq. ^ by estimating the correction 
SEk to the pure Dirac dispersion E^ = iej. coming from 
the second term in the angular brackets in the limit of 
small |k|. Expanding the square roots for M£ <C e k -C fi k 
this correction reads 



Mi 



(10) 



and is, therefore, negligibly small for the momenta close 
to the r point. The weak mixing of the two surface states 
is an artifact of the 2D model in the sense that in a real 
STI slab this effect will be suppressed exponentially in 
the slab thickness d. This discrepancy can be in principle 
removed by considering a different form of Mk that would 
more strongly vanish in the neighborhood of the T point. 
Such a form is easy to write in the momentum space but 
would necessarily imply longer range inter-surface hop- 
ping in the real space. Since our objective is to propose 
and test a simple model we shall not pursue this issue 
here, although we point out some consequences of this 
weak mixing when warranted. 

A simple but instructive test of the model can be con- 
ducted by adding to Hamiltonian ^ a term 



6H V = diag(y, V, -V, -V). 



(11) 



Physically this corresponds to introducing an opposite 
bias for the two surfaces which can be achieved e.g. by 
placing the sample inside a capacitor. Such a bias will 
play an important role in our considerations of the exci- 
ton condensate later in the paper. The spectrum of the 
modified Hamiltonian is displayed in Fig.[T|D. We observe 
that while the high-energy bands remain unchanged the 
Dirac points at low energies are simply moved up and 
down in energy, exactly as one would expect for the de- 
grees of freedom localized in a single layer. 

Another interesting modification consists of a replace- 
ment hk — > /ik 4- ms z . This represents a magnetic per- 
turbation with a moment along the z direction (i.e. per- 
pendicular to the surfaces) . Fig. [TJ; shows the spectrum 
that results from applying such a perturbation to one of 
the two surfaces. We observe that, as expected, one of 
the Dirac cones acquires a gap while the other remains 
gapless. 



B. Model II 

In some situations it will be desirable to simplify the 
surface Hamiltonian one step further. Specifically, if we 
are only interested in the low-energy surface degrees of 
freedom we may want to integrate out the four gapped 
bands appearing in the spectrum of Hi. A quick reflec- 
tion suggests that a very simple Hamiltonian of the form 



H n = 



M k 



will have the correct qualitative properties, if we take the 
two diagonal blocks to describe the two surfaces. The 
spectrum of this T-invariant Hamiltonian indeed has a 
single Dirac point at T associated with each surface. The 
low-energy properties are very similar to Hi and are il- 
lustrated in the bottom row of Fig. [I] 

These considerations suggest that the heuristically 
constructed lattice Hamiltonians ^ and (12 1 indeed can 
be used to describe the physics of a pair of surfaces in a 
slab of a STI. As mentioned previously, such a description 
offers a range of advantages in modeling such surfaces, es- 
pecially when quantitatively accurate predictions are de- 
sired. In the following Section we conduct further tests 
of these model Hamiltonians by comparing their spec- 
tral functions to the one derived from the exact surface 
propagator for the full 3D lattice model. 



III. EXACT SURFACE PROPAGATOR 

We now return to the idea of explicitly deriving the sur- 
face theory of a STI by integrating out the bulk degrees 
of freedom in a 3D lattice model. We start by writing 
down the action So for the fermionic degrees of freedom 
in a STI slab. At finite inverse temperature f3 — 1/fcgT 
this action reads 

rP _ 

So = / dT{V{d T + H S )V + X(9r + H b ) X 



+ ^Tt x + xT*}, (13) 

where r is the imaginary time, "J and \ are Grassman 
fields representing the surface and bulk fermionic degrees 
of freedom, respectively. H s is the Hamiltonian of the 
two surface layers, H^ describes all the bulk layers and T 
represents the terms that couple surfaces to the bulk. In 
the tight-binding description that we pursue in this work 
H s , Hb and T should be thought of simply as matrices in 
spin, orbital and lattice spaces; e.g. H s will be given by 
Eq. ([3]) for each surface. We are interested in the effective 
surface action defined as 



e -SM*,*] = V[x >X ]e 



-So[*,*;x,x] 



(14) 



(12) 



Since we are neglecting interactions at this stage the in- 
tegral is elementary and we obtain 

Scff= / dT^idr+H^-T^dr+Hb)- 1 ^} . 

Jo 

(15) 

For time-independent Hamiltonians it is useful to pass 
to the Matsubara representation by defining ^(r) = 
ft- 1 J2 n e~ iuj " T ^{iLo n ) with u) n = (2n 4- 1)tt//3 and n in- 
teger. The effective action then takes the form 

S ctf = ij^iw*) [Qj\iu n ) - Tig b (iuj n )T] *(iw„), 

n 

(16) 
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where Q a {iuj n ) = -{iw n -H s 
agator and Gb{iu n ) = —(itu. 
propagator. Eq. (161 informs us that the surface degrees 



1 is the bare surface prop- 
- Hi,)" 1 is the bare bulk 



of freedom are described by the effective propagator GeS 
of the form 



Ges{iv n ) = [G s 1 (iu n ) ~ T*Gb(Wn)T] 



(17) 



For a given 3D lattice model the effective surface the- 
ory can thus be obtained by performing the requisite ma- 
trix multiplications and inversions indicated in Eq. (17 1 



While this is certainly possible to accomplish in princi- 
ple, in reality one encounters at least two problems when 
implementing this procedure. First, for a system with 
many bulk layers Hi, and T will be large matrices (for 
the simplest model working in momentum space for each 
layer Hf, is a AN x 47V matrix where N is the number of 
bulk layers). Such matrices can only be inverted using 
numerical techniques. Second, G e s(iw n ) when analyti- 
cally continued to real frequencies, will not have a simple 
form — (oj — iJ c ff) _1 and therefore, except at low energies, 
will not yield a simple description in terms of a hermitian 
Hamiltonian. It is easy to see that for u; > A, the bulk 
bandgap, Gcff will have a complex- valued self energy, ex- 
pressing the physical fact that the surface electrons can 
decay into the bulk states which have been integrated 
out. 



A. Slab geometry 

In the following we shall perform a numerical evalu- 
ation of Geff hi a system with a finite number of bulk 
layers N and we also derive an analytic solution for G e B 
in the limit N —¥ oo and for special choice of model pa- 
rameters. We use these results to further validate our 
proposed effective Hamiltonian for the TI surfaces. In 
these calculations we use th e effe ctive model for TIs in 
the Bi2Se3 family of materialman regularized on the sim- 
ple cubic lattice. For a slab that is infinite in the (x, y) 
plane and consists of N + 2 quintuple layers (N bulk lay- 
ers and 2 surface layers) the relevant Hamiltonian can be 
written as a 4(iV + 2) x 4(N + 2) matrix, 



H 



/H R 

#t H R 

R} Ho R 

R) Ho 



\ 










\ 

Here Hq, defined in Eq. [3j describes each layer and 



H Q R 
i? f Ho) 



(18) 



R = - 



t z -X 2 
t z + \ z 



(19) 



is the interlayer coupling, with t z and X z representing the 
two types of interlayer hopping amplitudes permitted by 




FIG. 2: Spectral functions A(k, uj) for the surface states of a 
model TI as calculated from model I for a) orbital 1 and b) or- 
bital 2. Panels c) and d) show the exact spectral function for 
the surface layer calculated from Eq. (171 for a system with 
N = 46 bulk layers and again for orbital 1 and 2 respectively. 
For each orbital both spin projections are included. For the 
bottom row we use t — t z — 0.5, A z = 1.0 and e = 2.0 corre- 
sponding to bulk STI phase with Z2 index (1;000). Spectral 
functions are represented as density plots with light colors 
showing regions of high spectral density. A broadening pa- 
rameter S = 0.02 was used when producing these plots. 



the time-reversal and crystal symmetries. We note that 
R is proportional to a unit matrix in the spin space and is 
therefore a 4 x 4 matrix in the combined orbital and spin 
spaces. Here, and hereafter we treat the spin degrees of 
freedom implicitly. 

The surface and bulk Hamiltonians are now easily iden- 
tifiable from Eq. (18 1. We obtain H s = diag(ifo) Ho), 
T — (i?, 0,0,..., 0) and Hi, contains the remaining 
AN x 47V matrix. 



With these ingredients it is now straightforward to nu- 
merically evaluate GcS from Eq. ( 17) when the number of 



layers N is not too large. Fig. [2]shows the relevant spec- 
tral function A e ff(k, uj) = — ir Im(?efr(k, u) + iS) evalu- 
ated for N = 46 and a set of realistic model parameters 
indicated in the caption. The spectral function shows the 
characteristic gapless surface mode with the Dirac cone 
centered at the V point of the Brillouin zone and a con- 
tinuum of states above the bulk bandgap. As mentioned 
above this continuum reflects the physical fact that an 
electron injected into the surface at the energy above the 
gap has a short lifetime due to the hybridization with 
the bulk states. For comparison we also show the spec- 
tral function A(]&,uj) associated with one of the surfaces 
as calculated from our simple two-layer model I given in 
Eq. ([7]). The latter shows the same low-energy surface 
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state as the exact spectral function but of course cannot 
reproduce the continuum of levels that exists above the 
gap. Nevertheless we observe that the gapped bands that 
follow from Eq. ^ actually mimic the spectral form of 
the bulk states as well as can be expected from such a 
simple model. We also note that the low-energy exci- 
tations in both model I and the exact spectral function 
derive from a single orbital. We have also verified that 
the electron states exhibit the correct spin-momentum 
locking expected of the surface states in a strong topo- 
logical insulator. 



B. Semi-infinite slab 

As an interesting application of the formalism intro- 
duced above we now show how a simple analytic formula 
can be derived for Q c ff in the limit of semi-infinite slab 
N — > oo. This exact expression will be employed in our 
subsequent calculations of magnetic and exciton instabil- 
ities and can be useful potentially in other calculations. 
In this case we are concerned with a single surface and 
thus we take H s = H . To obtain Q c g in this situation 
we note that removing a single layer from the surface 
of a semi-infinite slab will necessarily leave the surface 
propagator unchanged. Since the electron hoping is only 
allowed between the neighboring layers it is permissible 
to describe the layer just below the surface by the effec- 
tive propagator Q c jf. This leads to the conclusion that 
Qb in Eq. (17) can be replaced by G c g. Thus, for semi- 



infinite slab we have a self-consistent equation 
Geg(iu n ) = [Go l {i^n) - T^Q eS (iu} n )T] 1 , 
where Qo(uo n ) = ~{iu n - -Ho) -1 - 



In its general form Eq. ( 20 ) involves solving for a 4 x 
4 matrix G e fi for each value of the momentum k and 
One can simplify this problem by first 



frequency uj„ 

noting that of all terms in Eq. ( 20 1 only has non-trivial 



dependence on electron spin {Mu and T are proportional 
to the identity matrix in the spin space). For each value 
of k one can thus perform a rotation in spin space — > 
= U^h^U^ with U = exp (— ijifi ■ s) and m the unit 
vector in the direction of (sin k x , sin k y , 0). The rotated 
Hamiltonian is diagonal in the spin space, 



Sz^k 



(21) 



In this 'helical' basis the problem for up and down spins 
decouples and we obtain two independent copies of Eq. 
(20) for the two spin projections at each k. The resulting 
equation for the 2x2 matrix Q c ff can now be solved 
analytically (it is, basically, a quadratic equation for a 2x 
2 matrix) . Unfortunately, for a generic set of parameters 
the explicit solution is quite complicated and does not 
lend itself to a convenient analysis. 

For this reason we focus in the following on a special 
case defined by the condition t z — X z (but all remain- 
ing parameters arbitrary). As one can see by inspecting 




M xx 



c) 





M XX 



FIG. 3: Exact spectral functions A^{k, u)) for the surface 
states of a semi-infinite TI slab as calculated from Eq. ( |30[ ). 
In all panels t z = t = X z — 1.0 and 8 — 0.02. Panel a) 
shows STI in (1;000) phase with e — 4.0, panel b) shows 
a critical point at e = 2.0, panel c) displays a weak TI in 
(0;111) phase for e = 1.0 and panel d) shows a (1;111) STI 
with e = —4.0. Note that in all plots spectral functions for a 
single spin helicity are shown. Those of the opposite helicity 
can be obtained by replacing oj — > —ui in the existing plots. 



Eq. ( 19 1 at this special point a further decoupling in the 



(20) orbital space occurs. Namely, if we define an effective sur- 
face propagator Q^g for each of the two orbitals o = 1, 2 
in the surface layer then it is easy to see that the matrix 
Eq. ( 20 ) breaks down to two coupled equations for scalar 



propagators y^ s 



0(2) 



(22) 
(23) 



Here H± = {—iuj n ± et), Q = —2X Z and M is defined in 
Eq. (pb. These equations now have a simple solution, 



c (i) 



2Q 2 # 4 



(H+H- + Q 2 



M 



f2\ 



(24) 



±y/(H+H- +Q 2 - M 2 ) 2 - 4Q 2 H+H- 



(2) 

and similar for Q cS . We shall see below that only the 
solution with the minus sign has the correct analytical 
structure and thus represents the physical surface propa- 
gator. Also, it is important to remember that Q^g repre- 
sents an effective surface propagator for orbital 1 with all 
other degrees of freedom (including orbital 2) integrated 

(2) 

out, and conversely for Q cff . 
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We now study the spectral function (k, w) associ- 
ated with the propagator (24 1. It is particularly interest- 



ing to see how an odd number of surface states emerges 
from this simple propagator. Upon analytical continua- 
tion iuj n — > uj + iS we see that the denominator has poles 
when lu — ek. At low energies the only contribution to 
the spectral function comes from these poles (because 
the square root in the numerator is always real) . Thus to 
study the surface states we can evaluate the propagator 
at the poles, obtaining 



Q 2 - Mi 



-2Q 2 (w 



\Q 2 

- iS - 



(25) 



The corresponding spectral function can be written, (at 
low energies) as 



Q -Mi 



6{lu 




when 
when 



Ml 



< 



(26) 



In order to have a non-zero spectral weight near the T- 
point we require that M^ =0 < Q 2 , or 



|e - 4t| < 2t z 



(27) 



In addition, to obtain a single Dirac point, we require 
that the spectral weight is zero at (0, 7r) and (7r, 7t), which 
translates to 



posed models, we now study the simple case when short- 
range components of the intra-layer and inter-layer repul- 
sive Coulomb interactions are included. We model these 
with the following Hubbard-like interaction terms, 



H lnt = U x 



E 



n l,x n 2,yi 



(31) 



x, V 



where x, y label the position, spin (and orbital for model 
I) indices. The interaction parameters for different or- 
bitals will generically be different but for the sake of 
simplicity we limit ourselves to the above two-parameter 
model. Furthermore, we shall only consider on-site inter- 
actions and will hereafter drop the position index. As al- 
ready mentioned, the leading instability associated with 
each interaction is considered leading to exciton conden- 
sation for the inter-layer interaction, and the ferromag- 
netic instability for the intra-layer interaction. We use 
a simple mean-field decoupling to diagonalize the lat- 
tice models exactly and the resulting phase diagram is 
obtained in the U1-U2 plane. The calculation carries 
through essentially unchanged from the equivalent calcu- 
lation for the continuous Dirac Hamiltonian^ with which 
we compare our results. 



A. Exciton instability 



> 2t 2 



|e + 4t| > 2t z 



(28) 



These conditions are easiest to analyze in the special case 
when t z = t. Then, a single Dirac point emerges when 



2t < e < 6t. 



(29) 



We note that this is exactly the condition for the model 
to be in the (1;000) topological phased i.e. STI with 
a surface Dirac fermion at k = 0. More generally, it is 
easy to see that by tuning e and t z one can access all 16 
topological phases that can be present in a 7~-invariant 
band insulator. 

In Fig. bl we plot the full spectral function (k, ui) 
for a selection of parameters representing some of the 
topological phases. When deriving the spectral function 
it is essential to pick the correct complex branch of the 
square root in Eq. (24). The retarded propagator that 



yields a positive definite spectral function and decays as 
~ 1/w at large frequencies has the following form, 



t (k,w) 



sgn(w)v/- 



g 2 +^Q 2 [{u + i5f 



y eff,retV 

with g = (10 + iS) 2 



-2Q 2 (uj + i8 

Q 2 



Ml 



(30) 



When including the inter-layer Coulomb interaction, 
electrons and holes on opposite surfaces can form exci- 
tons and condense. The interaction can be decoupled 
with a matrix-valued order parameter T = ^(^l^l) 
with the surface spinor in surface £, having 4 and 
2 components for model I and II respectively. Using 
the 8 or 4-component spinor for the combined surfaces 
\& = ($1, ^2) T , we write the expectation value taken 
with respect to the mean-field Hamiltonian as 



H I/II - rtl/ll 



r 
rt 



1 



* + — tr^r). (32) 



U, 



For chemical potential /j, close to zero, an order param- 
eter that opens a gap in the spectrum lowers the overall 
kinetic energy and will be favored. The maximal gap 
is obtained by considering an order parameter that anti- 
commutes with the kinetic part of the Hamiltonian. This 
form is diagonal in spin space and also preserves time- 
reversal symmetry. For model I we find 



Tr = 



2iAl\ 
-2iAl J 



(33) 



and for model II 



IV. INTERACTION EFFECTS 



zAl, 



(34) 



So far we have restricted ourselves to non-interacting 
Hamiltonians. As a further demonstration of the pro- 



with A assumed real and 1 being the identity matrix in 
spin space. The mean-field Hamiltonians can be diago- 
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nalized exactly to find the following spectra 



E 



CD 

k, s, a, b 



-/i + a 



V 2 + Ml + 2Vil + 4A 



2bJ(se k V -Vl) 2 + M 2 {V 2 + Q 



1/2 



and 



E, 



k, s, a 



(V + se k y+M 2 + A 2 , 



,(35) 



(36) 



with s,a,b = ±1. The relation between A and U 2 can 
then be obtained from the gap equations which result 
from minimizing the ground state energy 

4 J) =E' ^lab + ^N, (37) 



and 



Ei^=y' E^ +2—N. 



(38) 



where indicates a sum over occupied states and N is 
the number of sites of the 2D lattice. The resulting gap 
equations are 

U 2 \ V A 

A = ~8N^Ks,a, b ^m 77 (39) 

^k, s,a,b ^ V 



and 



A = -^V' A 

' " T ^k,s,a 



4 AT- 



p( J/ ) -a-,,' 



(40) 



These gap equations can be solved numerically for a 
finite, but large enough system (here a square lattice of 
1000 by 1000 sites) with a fixed number of electrons, or 
fixed filling fraction, as opposed to a fixed chemical po- 
tential. The number of electrons in a calculation is given 
by 8L 2 n and AL 2 n for model I and II respectively, with 
n G [0, 1] being the filling fraction and L 2 the number of 
sites. We present results for the half-filled case n = 0.5 
and for a slightly higher filling with n — 0.55 for model 
I and n = 0.6 for model II, such that we are compar- 
ing both models for the same number of electrons in the 
Dirac bands. Figures [4] (a,c) show half the size of the 
gap (the actual gap being 2A) associated with the EC 
instability of each model. 

Figures [4] (a,c) help illustrate an important point re- 
garding the exciton condensate in this system. For V =/= 
and zero chemical potential the latter formally occurs 
at an infinitesimal coupling U 2 because the relevant gap 
equations are of the BCS form, familiar from the theory 
of superconductivity. (At V = the density of states at 
the Fermi level Np vanishes so the formation of EC re- 
quires a finite critical coupling U 2c .) Nevertheless we ob- 
serve from Figs.[4](a,c) that for U 2 < U 2c the exciton gap 
rapidly vanishes even for a relatively large bias V = 0.5, 



consistent with the BCS exponential form A ~ e 



-1/U 2 N F 



Thus, as a practical matter, it would appear that biasing 
the system does not significantly improve the chances for 
experimentally observing the exciton condensate, unless 
the interaction strength can be tuned to the vicinity of 
U 2c . We also emphasize that the precise numerical value 
of U 2c obtained in this work should not be taken as a reli- 
able guide for experiments because of our reliance on the 
simple MF theory and our neglect of the long-range part 
of the Coulomb repulsion. Accurate determination of the 
critical coupling and the relevant critical temperature T c 
is a problem of considerable complexity with predictions 
for the latter ranging from mK to room temperatureP^H 
Although the MF calculation employed in this study can- 
not reliably predict the absolute magnitude of the exci- 
ton gap or T c , we expect that it is capable of indicating 
trends, as discussed above, and help in understanding 
the competition with other ordered phases, such as the 
surface magnetism discussed below. 



B. Ferromagnetic instability 

We mentioned that a magnetic moment on one surface 
along the z direction will open a gap in the Dirac cone as- 
sociated with this surface by breaking time-reversal sym- 
metry. Such a gap can also appear spontaneously due to 
intra-layer Coulomb interaction. We again seek a matrix 
order parameter that will decouple this interaction. We 
consider here the specific case where the Coulomb inter- 
action has the same strength on both surfaces such that 
both surfaces are in the same phase and the magnetic 
gap will open simultaneously in all Dirac cones. The or- 
der parameter will now involve only terms diagonal in 
the surface inde. 

We start with the more complicated case of Hamilto- 
nian I. For simplicity we will consider the Coulomb in- 
teraction to be of the same strength for all combinations 
of quantum numbers (layer, orbital and spin). Using n 
as the number operator, 



Hf 1 = Ui 



E 



s, t 



,(41) 



with the spin degrees of freedom s, t € {t, I}, 1 and 2 la- 
beling the orbitals and I being the surface index as before. 
We rewrite the mean-field version of this Hamiltonian by 
defining the average magnetization within an orbital o as 
the imbalance between spin up and down electrons 



mi, = y[(n«,(i,t) - (nt,o,i)], 
and the average number of electrons in an orbital 



Ui 



(42) 



(43) 
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FIG. 4: Gap associated with the EC and magnetic instability as a function of the strength of the short-range Coulomb 
interaction for Hamiltonian I (left column) for filling fractions n = 0.5, 0.55 and Hamiltonian II (right column) for filling 
fractions n = 0.5, 0.6. Panels (a) and (c) show half the EC gap as a function of the inter-layer interaction TJi and panels (b) 
and (d) show half the magnetic gap as a function of intra-layer interaction U\. Results obtained by solving the mean-field gap 
equations p9l, Eol, |48| and |52| for a finite system of 1000 by 1000 sites. 



The resulting mean-field Hamiltonian is 



T-Lj 



FM 



Hi 



1 1 



h, 2) 









D t t - to. 2 s : 



with 



Ce = [ng + m 2] and De = [ni + rii i 



(44) 



(45) 



where we defined fig — 1 +fit t 2 ■ The part of the Hamil- 
tonian proportional to fig is just a shift in the chemical 
potential at each surface and can be dropped here. At 
exactly half-filling and zero bias, ng. 1 = ng. 2 such that 
Ci, Dg and terms proportional to fif x and nj 2 can all 
be similarly ignored. Since we assumed the same interac- 
tion in both bands, we get a maximum gap opening when 
|ttz£. 1 1 = I to^.2| = to and the order parameter is simply 



T, 









-s. 



(46) 



This form can also be obtained by requiring that the or- 
der parameter anticommutes with Hj for V = 0. This 
holds for any filling fraction as the anti-commutation re- 
lation does not depend on this parameter, but to loses 
its meaning as the average magnetization for each orbital 
away from half-filling. For non-zero biases, no order pa- 
rameter that anticommutes with the Hamiltonian can be 



found and it is therefore possible that a different form 
will be preferred. We shall nonetheless use this same 
form even for finite bias for the sake of simplicity. As 
noted before, Eq. (44) can now be rewritten by replacing 
hk hu+msz with the dispersion relation remaining the 



same as Eq. (35 ), but with e| = 4(sin k x + sin k y ) + m 2 . 



The ground state energy and the resulting gap equation 
are 



- V' E {1) + 2 



and 



'an 



- y' 



m 



k, s, a, b tp(I) 1 

fi k. s. a, b M 



(47) 



(48) 



bV(V - sQl /e k ) 



(se k V-nl) 2 + M£(V 2 + Sll)] 



Model II with its single orbital is much simpler. Fol- 
lowing the same procedure as for model I we find that 
the interaction term 



-2,/int 



(49) 



can be decoupled for zero bias with the order parameter 

T /7 = tos z , (50) 

which anticommutes with Hn- m now preserves its 
meaning as the average magnetization for any filling frac- 
tion at zero bias. We chose again to use the same order 
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FIG. 5: Phase diagram as a function of the intra-layer (Ui) and inter-layer (U2) short-range Coulomb interaction strength for 
Hamiltonian I (left column) and Hamiltonian II (right column). Main panels show true main region of interest while the insets 
sketch the full phase diagram. Panel (a) and (c) show the phase diagram for half-filling. Panel (b) shows model I at re = 0.55 
and panel (d) shows the phase diagram of model II at n = 0.6. The filling fraction of (b) and ( d) a re cho sen to have the same 
number of electron in the Dirac bands. Results obtained by solving the gap equations (39 1, (40 1, (481 and (521 for a finite 
system of 1000 by 1000 sites. 



parameter even at finite bias for simplicity even if a differ- 
ent form could be favored. Using the same substitution 
for ek, the spectrum is still given by Eq. (361 while the 
ground state energy is 



9 £-^l\r. S .„ 1 



k, s, a 



m »T 

H N, 

Ui ' 



(51) 



and the gap equation 



m 



■ C/ k, s, a 



-(l + sV/e k ). (52) 



The magnetic gap equations obtained are solved nu- 
merically in the same fashion as for the EC instability 
and panels (b) and (d) of Fig. [4] show half the size of the 
gap (the actual gap being 2m). 



C. Phase diagram 

We can now draw the phase diagram of a STI as a func- 
tion of the strength of both the intra-layer and inter-layer 
interaction. The phase boundary between the magnetic 
phase (MP) and the EC is calculated by evaluating the 
ground state energy and the coupling strength as a func- 
tion of the gap simultaneously to obtain the functions 
Ui(E g \A=o) an d U2(E g \ m= o) and plotting the resulting 
parametric curve in the U1-U2 plane. Fig. [5]shows the 
resulting phase diagram in detail while Fig. [6] compares 
model I and II to the continuous Dirac Hamiltonian ([I]). 



It should be noted that all models show a critical intra- 
layer coupling below which we are in a non-magnetic 
phase. For the EC however such a critical inter-layer 
coupling should only exist for zero bias and otherwise 
vanish. Fig. [2] and [5] show that the lattice models in- 
stead always have a critical inter-layer coupling (shown 
as dashed lines) even at non-zero bias. Below both criti- 
cal couplings we are in an uncondensed, or semi-metallic 
(SM) phase. This is an artifact of the lattice model due to 
a very small gap opening at /i = for non-zero bias. The 
size of this gap increases slower than ~ O.IV^ 3 for model 
I and ~ 0.25V^ 2 for model II in the parameter range con- 
sidered here. This only becomes relevant when studying 
the weakly interacting case where the continuous model 
predicts an exponentially small gap, and therefore diffi- 
cult to observe experimentally. In the region of interest 
for experiments and possible devices, this small gap can 
be ignored safely compared to the EC gap. 

Fig. [6] with its comparison of model I and II to the 
continuous Dirac Hamiltonian (I]) raises a number of in- 
teresting points. Although the low-energy physics of the 
surface states of a STI is generally thought to be well de- 
scribed by this effective continuous model, a comparison 
to our own effective lattice model suggests that it will be 
insufficient for a quantitative study of real TI materials. 
Not only does depend on an arbitrary cutoff, the exact 
position of the phase boundary differs significantly at in- 
termediate couplings. This means that estimates derived 
from the Dirac Hamiltonian for the critical temperature 
of the EC, for example, are unlikely to be very accurate. 
A more ambitious study with any hope of achieving a 
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FIG. 6: Comparison of the phase diagram for the continuous Dirac Hamiltonian |TJ and the lattice Hamiltonian I in panel 
(a) and Hamiltonian II in panel (b). This phase diagram is for the half-filling case. The momentum cutoff A used for the 
continuous model (shown in the insets) is tuned to reproduce the critical intra- layer coupling Ui c obtained from the lattice 
models I and II. Results obtained by solving the gap equations (391, (40 1, (48 1 and (52 l for a finite system of 1000 by 1000 
sites. 



quantitative estimate will require a model such as those 
proposed herein, that accounts properly for the underly- 
ing lattice structure of the STI and its more complicated 
band structure. 

Despite its limited applicability, the mean-field calcu- 
lation above also allows us some general observations of 
use for the experimental observation of EC. Fig. [5] sug- 
gests for example that an appropriate system for the ob- 
servation of a EC requires an inter-layer coupling at least 
twice as strong as the intra-layer coupling, otherwise the 
magnetic phase takes over as the most favorable state 
and the topological properties of the surface states are 
lost due to the breaking of time-reversal symmetry. 



D. Comparison to 3D results 

To ascertain the accuracy of the proposed models in 
terms of predicting the correct phase diagrams we now 
compare the results obtained above to a calculation based 
on the exact surface propagator derived in Sec. III. We 
remark that employing the exact surface propagator is 
equivalent to performing the calculation within the full 
3D lattice model. 

To this end it is useful to reformulate the mean field 
theory in terms of the path integral of Sec. III. Thus, we 
write the action for the 3D interacting TI as 



S* mt = So 



(53) 



where So is the noninteracting action d efine d by Eq. ( 13 ) 
and H mt represents the Hubbard term (31 1. For simplic- 



ity we consider only the on-site interaction U\ and neglect 
interactions between different orbitals. On every site the 



interaction term can be decoupled using the Hubbard- 
Stratonovich transformation based on the identity 



e 2 



where C is a constant and m represents an auxiliary real 
field. Noting further that (n-j- — n^) 2 = iVf + nj, — 2n^n^ 
one can express the Hubbard term in terms of the total 
charge and the total magnetization on each site, thus 
rendering the action ( 53 1 formally bilinear in Fermionic 
fields, 



S[m] = S 



rP 

o ^ 



2U X 



+ m a (r)(n a f - n aJ .) 



(55) 

where a represents the spatial coordinates and orbital 
indices. The (n-j- + n-j.) term from the decoupled H lnt has 
been absorbed into the chemical potential term. 

The mean field approximation consists of neglecting 
the temporal (quantum) fluctuations in m a (r), which can 
then be interpreted as on-site magnetization. If we fur- 
thermore assume a homogeneous solution, i.e. m a = m 
for all sites and orbitals, we can pass to the Fourier space 
and the magnetization term is then seen to modify the 
non- interacting action So by replacing — > + ms z . 
The equilibrium magnetization then follows from mini- 
mizing the effective action S e[m] with respect to m. The 
former is obtained by integrating out the Fermi fields in 
S[m], 



-SeffM 



p[x,x;*,*]e 



-S[m] 



(56) 



We proceed in two steps. First we integrate out the bulk 
degrees of freedom and focus on the surfaces. The pro- 
cedure is identical to the one outlined in Sec. III. For a 
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FIG. 7: Comparison of model I, II with the exact result for the 
dependence of magnetization m on the interaction parameter 
Ui. The parameters are t z = t = \ z — 1.0, e = 4.0 which 
corresponds to the bulk STI in (1;000) phase and V = fj. = 0. 



semi-infinite slab geometry it leads to the effective surface 
action given in terms of the surface propagator de- 



fined in Eq. (|24|) with £k 
where s — 



sW4(sin k x 



sin 2 k„ 



now labels the two spin projections ex- 
pressed in the helical basis. The effective MF action is 
thus given by 



-Sett [m] 



(57) 



exp 



y ' n,k,s 



0eff ( iw njk, S 



-* s (icj„,k) 



The integral over the remaining Fermi fields is easily eval- 
uated and leads to 



S oS [m]=P — N 



E m [Sis ( iuJ n,K s) 



(58) 



n,k,s 



Before tackling the minimization it useful to first per- 
form the summation over the spin index which yields a 
more transparent expression 



S cS [m]=f3 — N+ y > 'In 



n,k 



D 



(59) 



quite good, especially for model II whose critical coupling 
U([ falls within 2% of the exact value U c = 3.39 for the 
model parameters specified in the caption of Fig. [7] We 
have considered other values of model parameters and 
found similar level of agreement, with the model I and 
II critical couplings lying typically within 10-15% of the 
exact value. In all cases we found that model I and II 
critical couplings exceed the exact value. 



V. CONCLUDING REMARKS 

The key advance presented in this work consists of the 
formulation of two simple 2D lattice models that faith- 
fully describe the low-energy electronic states associated 
with a pair of surfaces of a strong topological insulator. 
We have demonstrated that these models produce elec- 
tron states with properties characteristic of true STI sur- 
face states, in terms of their spectral properties as well 
as their susceptibility to the formation of ordered phases 
in the presence of interactions. Employing these surface- 
only models will confer a significant advantage on any 
numerical calculation by obviating the need to explic- 
itly consider the (uninteresting) bulk degrees of freedom. 
This includes computations that seek to address the ef- 
fects of interactions (both electron-electron and electron- 
phonon) as well as those of various types of disorder or 
other perturbations, such as magnetism or superconduc- 
tivity, acting on the STI surface. 

Although we have focused here on the example of two 
parallel surfaces one can imagine applying the same strat- 
egy to other geometries and situations. One example that 
we have considered (but not discussed here in detail) is 
the surface of a cylindrical STI wire which can be mod- 
eled by 2D lattice Hamiltonians similar to those defining 
our model I or II with an extra hopping allowed between 
radially opposed sites on the cylinder. Also, one can 
straightforwardly apply this strategy to 2D topological 
insulators (as well as 2D Chern insulators) and formu- 
late ID lattice models for the pair of edges. Finally, one 
can imagine constructing lattice models for 2D surfaces 
or ID edges for 3D or 2D topological superconductors 
using our techniques since their first-quantized, single- 
particle Hamiltonians are in many cases identical to those 
describing TIs. 



with D n = Q^) (icj n ,k)-ff_|_. Differentiating the result in 



Eq. ( 59 ) with respect to m leads to the exact gap equation 
(which we omit here for the sake of brevity). It is a 
simple matter to solve this gap equation numerically and 
compare the outcome to our previous results based on 
the surface-only models I and II. The dependence of m on 
U\ is shown in Fig. [7] We observe that the agreement is 
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